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We report the synchronization behavior in a one-dimensional chain of identical limit cycle oscilla- 
tors coupled to a mass-spring load via a force relation. We consider the effect of periodic parametric 
modulation on the final synchronization states of the system. Two types of external parametric ex- 
citations are investigated numerically: periodic modulation of the stiffness of the inertial oscillator 
and periodic excitation of the frequency of the self-oscillatory element. We show that the synchro- 
nization scenarios are ruled not only by the choice of parameters of the excitation force but depend 
on the initial collective state in the ensemble. We give detailed analysis of entrainment behavior 
for initially homogeneous and inhomogeneous states. Among other results, we describe a regime 
of partial synchronization. This regime is characterized by the frequency of collective oscillation 
being entrained to the stimulation frequency but different from the average individual oscillators 
frequency. 

PACS numbers: 05.45.Xt,05.45.Pq,02.30.Mv 
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INTRODUCTION 

In Nature many life significant processes are regulated 
by mechanisms of synchronization and entrainment. A 
properly chosen external stimulus can synchronize an 
intrinsic biological cycle. A great number of biologi- 
cal processes involve several oscillatory cycles that in- 
teract and contribute to an overall system temporal be- 
havior. To give a few examples, the human heart re- 
sponds sensitively to the forced oscillations modulated 
by a pacemaker with various frequencies fH , the respira- 
tory rhythm in humans and animals can be entrained to 
a mechanical ventilator phase Fibrillar flight mus- 
cles of some insects act as a mechanical resonator that is 
coupled to the insects thorax and wings Q . A rhythm of 
flight muscles contractions and wings vibration frequency 
can be turned by an experimental mechanical or optical 
stimulation [1, [a]- In all the examples provided so far 
a self-oscillatory rhythm is modified due to mechanical 
forces via the coupling to an external mechanical oscil- 
lator. A simple theoretical description in terms of in- 
teracting nonlinear oscillators would provide a universal 
approach for exploiting dynamics in this type of systems. 

For analysis of time evolution processes simplified 
models from nonlinear dynamics have been used [6|. In 
fact, in the framework of nonlinear models one can de- 
scribe the phenomena that occur in a relatively large class 
of biological and chemical systems 0, Q . In this paper we 
use the analogies found in nature and introduce a non- 
feedback model of entrainment in a system of coupled 
nonlinear oscillators. Our model might help to elucidate 
the basic underlying mechanisms that regulate synchro- 
nization of processes in a real biosystem. 

Although simplified low-dimensional models can cap- 
ture inherent mechanisms leading to synchronization, 
quite often it happens that spatial dynamics plays a cru- 



cial role in mediating synchronization fo'] in a given sys- 
tem. To unveil the whole variety of possible synchroniza- 
tion scenarios one has to introduce a spatial dimension 
in the mathematical model of a biosystem. 

Collective phenomena in spatially extended systems 
and their one-dimensional chain analogues have been a 
subject of a large body of investigations [^, [3, For 
instance, models of self-oscillatory systems are widely 
used to describe biochemical pattern formation processes 
in spatially extended systems [13, • It is well known 
that the complex spatiotemporal behavior in a system is 
largely determined by its dynamical instabilities. There 
are two basic approaches that are used for numerical and 
analytical studies of stability in spatially extended sys- 
tems subjected to external forcing. The first approach 
implements the mathematical framework provided by the 
complex Ginzburg-Landau equation (CGLE) [l^. In 
this case the description of an oscillatory medium is re- 
duced to the study of the amplitude equation with a 
periodic forcing. Different types of periodic forcing of 
the CGLE and a rich set of dynamical states and their 
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bifurcations have been studied in 
second approach uses a description of system spatiotem- 
poral dynamics in terms of a large populations of coupled 
oscillators. Due to a large number of degrees of freedom 
these systems can display a pano ply of dynamical behav- 
iors: from cluster formation 16..Tl7l.[l8l| to multistabilitv 
regimes [l^, [3| with the coexistence of distinct stable 
collective modes of oscillations. 

Most of the theoretical studies mentioned so far are de- 
voted to investigation of dynamics in globally or locally 
coupled ensembles of oscillators. The coupling term usu- 
ally depends on phases or displacements of oscillators. In 
this work we discuss a new type of coupling via a common 
force. This force explicitly depends on the displacements 
as well as on the velocities and the accelerations of indi- 
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vidual elements in a large ensemble. The idea to intro- 
duce such a coupling has been first proposed in Ref. [1^ 
The model of a large ensemble of self-oscillatory elements 
coupled via a common force to an inertial oscillator was 
proposed in [2^ as a phenomenological approach to study 
spontaneous oscillations in an active muscle fiber coupled 
to mechanical resonator 2l|. Each element in the en- 



semble is described by the Stuart-Landau equation with 
real coefficient. It is shown that the system undergoes 
a Hopf bifurcation near to a critical threshold defined 
in systems parameter space. Various regimes of collec- 
tive oscillations are reported to exist. In fact, the sys- 
tems parameter space is divided into several regions with 
mono- and multi-stable collective modes of oscillatory 
dynamics. Near to an instability threshold defined in 
the parameter space there exist completely synchronous, 
asynchronous and antisynchronous regimes of collective 
dynamics. In addition, regions with the coexistence of 
several modes are found. 

In this paper we introduce a modification of the model 
considered in Ref. [13]. We include an external time- 
dependent stimulation by adding small amplitude time- 
dependent terms in the parameters expression. We de- 
scribe the mechanism of entrainment of the frequency of 
system solutions to the stimulation frequency. The es- 
sential feature of the described non-feedback mechanism 
of entrainment is that the external stimulation enters in 
the model as a parametric modulation. Two types of 
parametric excitations are treated: the first type is the 
periodic modulation of the inertial oscillator stiffness; the 
second type is the peridic excitation of a natural fre- 
quency of a self-oscillatory element. We show that the 
synchronization scenarios are ruled not only by the choice 
of parameters of the excitation force but depend on the 
initial collective state in the ensemble. For both types of 
the periodic excitations entrainment behavior is studied 
for homogeneous states. We find sets of frequency-locked 
solutions that correspond to the resonant driving. We 
also study the effect of two types of periodic driving on 
stability of inhomogeneous states. In this work we ad- 
dress several questions: what is the influence of periodic 
parametric forcing on the stability of a particular collec- 
tive motion, and how does the features of synchronization 
compare for different types of parametric excitations in 
the system. 

This paper is organized as follows. In Section , we in- 
troduce the autonomous model and briefly recall some 
features of the complex phase dynamics from Next, 
in Section we provide the time-dependent model with 
the first type of parametric modulation and investigate 
the existence of stable periodic solutions for the case of 
A'^ = 1 oscillator in the chain. We explain the numeri- 
cal procedure that is used for the amplitude-frequency 
response calculations. Amplitude-frequency responses 
from numerical simulations and analytical derivation are 
compared. Next, we present, using numerical simula- 



tions, the effect of the periodic forcing on initially inho- 
mogeneous states behavior. We list a variety of regimes 
of inhomogeneous collective dynamics including partially 
synchronized states. The second type of parametric mod- 
ulation is introduced in Section . Similarly to the pre- 
vious section, we provide stability diagrams for the case 
of A'^ = 1 oscillatory unit in the chain coupled with an 
inertial element. We discuss a biological relevance of our 
findings and give general conclusions in Section . In the 
Appendix an analytical amplitude-frequency relation is 
derived using a first order quasiperiodic approach . 



MODEL EQUATIONS 

In this section we review the basic model [l^l of a chain 
of active mechanical oscillators coupled to a damped 
mass-spring oscillator. Each elementary unit in the chain 
can be modeled by the Stuart -Landau equation. Two 
variables are defined describing the amplitude r.^ and the 
phase of each ith element in a finite chain of N ele- 
ments: 



(t>i = 



(1) 



where B is the Landau coefficient. The set of N equations 
for i = 1^ . . . N describes N uncoupled oscillators. In this 
paper we consider mechanical elements all arranged in a 
one-dimensional chain and coupled via a common force 
due to a linear inertial oscillator attached to one end 
of the chain. In the model [20| a chain of N elements 
has one fixed boundary and is connected to a mass load 
in a way that the change of the total extension of -|- 
1 oscillators is zero. This corresponds to the no flux 
boundary condition. Consequently, this condition can be 
recast in the form of the following geometrical constraint 
imposed on oscillators motion: 



N+l 
1=1 



(2) 



Where Xi = Ti cos (f>i,i = 1, . . . N is the displacement of 
the ith active self-oscillatory element, xn+i is a displace- 
ment of the linear inertial oscillator. The external force 
due to the mass-spring load is determined by the equa- 
tion of motion: 



N 



F 



(3) 



where the parameters M, SIq and v are the mass, char- 
acteristic frequency of linear oscillator and the spring 
damping. The spring stiffness coefficient is k = Aff^Q. 
Upon including the mechanical force F from ^ Eqs.((T|) 
can be written in the complex form : 



{iu + e)z,-B\z,\^Zi + FS,- 



(4) 
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where Zi = riC^'^^ and ^ is the scahng parameter. The 
set of equations ©-dH) represents a system of active os- 
cillators that are coupled via a mean field F/^. Unlike 
in the Kuramoto model of phase oscillators with global 
coupling lj| here the collective behavior is governed 



by a mean field that depends on the velocities and the 
accelerations of elements in the chain. It was shown in 
Ref.[23| that oscillatory versus non-oscillatory motions 
exist in an unfolded space of the system parameters. 
Namely, oscillatory versus non-oscillatory behavior de- 
pends on the values of stiffness, mass of inertial load and 
the bifurcation parameter e. In fact the onset of the os- 
cillatory behavior happens for e < when the system 
undergoes a Hopf bifurcation at the critical value of lin- 
ear frequency fig = f^c- The critical frequency for which 
transition from a stable non- oscillatory to an oscillating 
behavior occurs is defined in 20] from the expression: 
ilc = {uj^/{l + NMv/C) + e{iy-e + 2^/{NM))y/^. While 
oscillatory states can be found for fio < for higher 
values of the linear frequency f2o > initial oscillations 
always decay to the stable trivial solution. 

A rich bifurcation diagram with various modes of col- 
lective oscillatory dynamics exists for e > near the crit- 
ical instability threshold [2^ . Several regions with homo- 
geneous and inhomogeneous collective modes as well as 
regions with coexistence of different dynamical regimes 
are found. Depending on the frequency and e one can 
observe different asymptotic states of system behavior: 
asynchronous, antisynchronous and synchronous states. 
For a low value of the stifness fio <C f^c the system is set 
into a synchronized oscillation state. The phase and the 
amplitude of self-oscillatory elements in a globally syn- 
chronized regime are identical. For large positive e and 
for r^o nearly close to the instability threshold value flc 
a set of asynchronous states is found. The asynchronous 
state can be characterized by an inhomogeneous distri- 
bution of phases and amplitudes of oscillators along the 
chain. Remarkably, that while the total chain extension 
vibrates at one frequency the value of the average os- 
cillators frequency is below it ^2d\. Localized regions of 
synchronized motion can be identified along the chain. 
For large value of the frequency flo and e > the system 
is in an antisynchronous regime. In this state neighboring 
oscillators move in an antiphase so that the total chain 
extension remains fixed. 



EFFECT OF PERIODIC MODULATION OF 
STIFFNESS 
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FIG. 1: Arnol'd tongues: stability diagram for the parameter 
space defined by the amplitude Ai and frequency ratio Qei^o^ ■ 
The boundaries that separate stable limit cycles from unstable 
and quasiperiodic solutions are shown. Frequency ratios of 
entrained limit cycle solutions inside each stability region are 
indicated. Parameters: e = —0.2, v = 0.2, ^ = 0.3125, B = 
1,LJ = 1, M = l,k = 0.36. Original unperturbed limit cycle 
frequency: loq — 0.7439. 
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FIG. 2: Examples of limit cycle solutions. Thick fine shows 
the original limit cycle for K{Xi,t) = fc in system Q-®. Two 
other limit cycle solutions are calculated for the system Q- 
^ with non-zero amplitude of driving Ai — 0.5 k (thin fine), 
and with Ai = fc (dash-dotted fine). External frequency fie = 
ujo ~ 0.7439 and the stiffness coefficient k = 0.36. 



In this section, we discuss the effect of a non-additive 
periodic excitation on lock-in solutions stability in the 
system of self-oscillators coupled to an inertial load. The 
excitation is taken in the form of periodic modulation of 
the inertial oscillator stiffness coefficient. First we define 
a set of equations of motion and then proceed to exam- 



ining the system stability domains. Wc study the am- 
plitude dependence for frequencies of excitation near to 
the one-to-one resonance. Let us introduce a parametric 
modulation in the stiffness k of the inertial oscillator ([3]) 
by using a periodic term q{Xi,t) = Aicosile^ with fre- 
quency Qe and amplitude Ai. Then the expression for 
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the force is: 



F : 



N 



+ K{\i,t)zi + vzi) 



(5) 



where K{\i,t) = fc + (7(Ai,i) is the time-dependent stiff- 
ness coefficient. In our treatment we include the case of 
a strong magnitude of driving by taking Xi/k < I. Upon 
differentiating both side of Eq. (|4|) and substituting the 
expression for Zi from ([5]) the differential equation for the 
force yields: 



F = - 



ML 

NM 



1 

N 



N 



^Re 



d(B\z,\^i 



dt 



^(e + v)zi - K{Xi,t)zt 



(6) 



Together with Eq. (|3|) we now obtain a complete set of 
differential equations for + 1 variables. The displace- 
ment xn+1 of the inertial oscillator is determined from 
the geometrical constraint ^ . We expect that the stabil- 
ity behavior in system ([4])- ([6]) is largely determined by a 
new type of coupling. As a consequence of the coupling 
the motion of every self-oscillatory element is affected 
by a common force. This force is explicitly dependent on 
time as well as on the oscillators displacement, velocity 
and acceleration. 



Entrainment for an initially synchronized regime 

Let us first discuss the entrainment behavior for ini- 
tially homogeneous solutions: all the oscillators in en- 
semble are moving in phase and with the same frequency 
of oscillations. The linear stability analysis of the time- 
independent system in Ref. [lO] provides a detailed phase 
diagram with the regions of stable globally synchronized 
solutions. Indeed, for the parameters space correspond- 
ing to e < and ilo < there exists a large set of 
stable synchronized solutions of Eqs. ([I])-® for Ai = 0. 
When the synchronized state is stable it suffices to re- 
duce systems description to the case of 1 oscillator 
coupled with inertial load. However, when the paramet- 
ric modulation is introduced in the model the stability 
problems for the system of = 1 and iV > 1 oscilla- 
tors in the chain are no longer equivalent. The applied 
periodic excitation can change the stability of the syn- 
chronous solution for > 1. Thefore, the system can be 
driven to a new asymptoticcally inhomogeneous state. In 
our work we do not analyze the stability of synchronized 
solutions for the entire parameter space. Instead, we con- 
sider entrainment behavior of the synchronized state for a 
particular choice of parameters: e — —0.2 and f2o = 0.6. 
We integrate numerically the time-dependent system (jlj)- 
© with A^ = 1 and A^ = 100 oscillators. It follows that 
the observed entrainment behavior for the low-(A^ = 1) 



and high-dimensional (A^ > 1) cases matches well for 
the choice of parameters given above. Therefore, in or- 
der to simplify our analytical derivations and to speed up 
numerical procedure we reduce dynamical description of 
the system to the low-dimensional case. In the following 
paragraphs we discuss entrainment behavior for A'^ = 1 
oscillator coupled with a mass load. For initial conditions 
we chose arbitrary small displacement xi = zi. Results 
of numerical simulation are given in Fig. ([T]). We show 
stability domains in the parameter space defined by the 
amplitude of driving Ai and the ratio between external 
frequency fig and the original limit cycle frequency wq. 
Here the stability regions with lock-in solutions are sep- 
arated from the regions of quasiperiodic and unstable 
solutions by defined boundaries. Inside each stability re- 
gion the lock-in solution has rational frequencies ratio 
m : n. Note that the width of entrainment regions short- 
ens as one proceeds from the 1 : 1 resonance to the lower 
frequency ratios. Several limit cycle solutions are pre- 
sented in Fig. ^ for the 1 : 1 resonance and for different 
amplitudes Ai. 



Amplitude— frequency response 

In this section we present numerical and analytical 
results on the amplitude-frequency dependence for the 
system (j4l)-([6]). We use numerical analysis of funda- 
mental frequencies [2^ to determine the amplitude of 
quasiperiodic and periodic solutions of the system (|4])- 
^ for the frequencies close to the one-to-one stimu- 
lation. For solutions x = Rezi and y = Imzi of (j4|)- 
(l6|) we obtain a quasiperiodic approximation near the 
vicinity of the original limit cycle solution of Eqs. 
^ for Ai — 0. By using the substitution of solutions 
x{t) and y{t) with their zero time-averaged equivalents: 
x{t) x{t) — X and y{t) — > y{t) — y we write down the 
quasiperiodic expansions for the new x{t) and y{t): 



(7) 



where is a set of time-independent frequencies. We as- 
sume that in the above expression the amplitudes Cq and 
dq do not dependent on time. In the above expressions 
the terms are collected according to decreasing order of 
magnitude. To begin our numerical procedure we intro- 
duce a set of new coefficients c„ and dn that are expressed 
from the following integrals [S^l : 



1 



5n = ^ / x{t)e-^"^-^'>'xm, (8a) 



dn = ^ / j/(t)e-™(-/^)*x(i)di, (8b) 



2T 

where x(t) = 1-1- cos nt/T is the Hanning window. The 
interval of time IOOtt <T< SOOtt is the integration time 
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(a) 




FIG. 3: Amplitude-frequency response surface is shown on 
the (ilejAi) plane, (a) Numerically calculated amplitude 
ci of the leading frequency component of the limit cycle in 
Eq. ([H near to the 1 : 1 resonance for the driven system 
(S])-®. Stability boundaries of the 1 : 1 Arnol'd tongue 
are shown on the plane. (b) Numerically esti- 

mated maximum of the limit cycle solution for (|4])-(|6]) with 
max|3;| — max2007r^t^3007r \x{t) — x\ for different amplitudes 
Ai and frequencies fie of external driving. Parameters are: 
e = -0.2, u = 0.2, k = 0.36, M = 1, ^ = 0.3125. 




FIG. 4: Amplitude-frequency response curves obtained at a 
fixed amplitude of A = Xi/k = [0.7,0.8,1] from Fig. Q. 
Analytical amplitude-frequency dependence r^^(f2e) is dis- 
played (small circles) (see appendix for derivation). Here 
Q,eLdQ^ — Q,Qq^{1 — A), where A = 0.0711"^ is the frequency 
mismatch. Parameters are given as in the Appendix. Numeri- 
cally evaluated stability boundaries are shown (crosses) for all 
three values of A. Parameters for the analytical amplitude- 
frequency dependence: e = —0.2,v = 0.2, M = 1,uj = l,f = 
0.3125, fio = 0.7861, Ai = Mflg 



for every systems trajectory. Next, two sequences of co- 
efficients Cq, and dq are produced by rearranging c„ and 
dn according to the decreasing order of magnitude. In 
the expansions ([7]) we only retain the first two terms 
that produce reasonably accurate approximation to the 
solutions x{t) and y{t). We calculate the first four largest 
coefficients of the series Cq and dq as follows: 

ci-i = max{ci,_i, C2,-2, ■ ■ •}, (9a) 
= max{Ji,_i,(i2,-2, • • ■}, (9b) 

Figure ([3^) illustrates the parametric dependence of 
the largest coefficient ci in the expansion ([5^) on the am- 
plitude Ai and the frequency Qe- Numerical responses 
are obtained by taking initial conditions for x{t) and 
y{t) to be randomly distributed in the vicinity of un- 
perturbed limit cycle for Eqs. ([I])-®. For comparison 
we also calculated the maximum of the numerical solu- 
tion over the same time interval T. Figure ^jp) shows 
numerical results for the maximum of the amplitude 
x{t): max|x| = maxT|x(i)|. In the plane (rJe,Ai) sta- 
bility boundaries for the 1 : 1 resonance are indicated. 
The results from two plots are consistent except for the 
regions near to the stability borders. The amplitude- 
frequency response is a surface that is smooth inside the 



borders of stability region and shows discontinuities out- 
side near to the borders of Arnold tongue. Apparently, 
the discontinuities are due to the emergence of instabil- 
ities. Because of the instabilities, the leading frequency 
components that are defined from the quasiperiodic ap- 
proximation ([7]) in this case might be distinctly different. 
Consequently, the numerical procedure produces a set of 
amplitudes which arc not smoothly dependent on the ex- 
ternal frequency fig- 

We compare numerically estimated amplitude- 
frequency responses with the results obtained by using 
quasiperiodic approximation (see Appendix). In Fig. ^ 
three numerically calculated amplitude-frequency curves 
are displayed for different values of the amplitude 
A — Xi/k. To compare these results with the analytical 
calculations, we display the amplitude-frequency curve 
rj^^(rie) obtained from the quasiperiodic approximation 
of hmit cycle solution {x{t),y{t)) (see Eq. (fT8|) in the 
Appendix). It is apparent, that for a given choice of 
parameters the curves from numerical and analytical 
calculations follow similar nonlinear response behavior 
inside the stability region (stability boundary is marked 
by crosses). 
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FIG. 5: Space-time diagrams of intensity of oscillations 
\zi\ for the ith element, (a) Asynchronous state for time- 
independent system (Ai = 0). (b) Asynchronous 
state for the time-dependent system Q-® (Ai — fc).The fre- 
quency of driving f2e = wo = 0.85. .(c) Synchronized behavior 
is obtained with the parametric forcing of an initially asyn- 
chronous state of system (I10|l - (|11| ) (A2 — 0.5a;, S7e — 0.85). 
Parameters for (a),(b) and (c): fio = 0.85, e — 0.5, k = 
Mnl,M ^l,uj = l,i/N = 0.3125, iV = 100. 



Entrainment for anti- and asynchronous regimes 

In this section we illustrate results of numerical analy- 
sis for the periodic forcing of initially spatially inhomoge- 
neous states. The existence and stability of these states 
for the autonomous system (Ai = 0) is controlled by the 
choice of a mass of inertial oscillator and parameter e. 
For the time-dependent system it is expected that under 
the influence of a strong driving these states might be- 
come unstable and converge to a new asymptotic state. 
In the limit of a large excitation it can produce signifi- 
cant changes on the regions of stability of existing spa- 
tiotemporal regimes. We demonstrate that, meanwhile 
some collective regimes do not survive new spatially ho- 
mogeneous states become asymptotically stable. Here 
we show which conditions on the frequency and the am- 
plitude of external driving have to be satisfied in order 
for the resulting spatial state to be stable and homoge- 
neous. In the literature studies of collective dynamics 
of chains of coupled oscillators are often reduced to their 
phase dynamics studies 0, [2^ . Our system is different 
from these cases because the amplitude of external force 
F strongly depends on the state of the ith oscillator, and 
therefore, on the collective dynamics in an ensemble. It is 
not sufficient to confine the investigation only to a phase 
dynamics study, like in the cases of weakly coupled oscil- 
lators [24*1. It should be pointed out, that in the limit of 
a strong excitation the effect on amplitude variation of 



FIG. 6: Frequency ratio curves for the antisynchronous (Anti) 
(asynchronous (A)) state of a chain of A'^ = 100 oscillators, 
(a) Ratio between external frequency fie and the average os- 
cillator frequency ujm = ^/N'^^ Ui is plotted versus external 
frequency for an initially antisynchronous state of the sys- 
tem (I31)-([Sl)-(b) Frequency curves are plotted for an initially 
asynchronous state of the system Q-®. Shown are two fre- 
quency ratio curves: Q.eLo^ versus He (squares) and Q,eijJ^ot 
versus fie (dots), where utot is the frequency of the total ex- 
tension Xtot of a chain of oscillators. Parameters: e = 0.1 
for antisynchronous state, (e = 0.5) for asynchronous state, 
Vlo = 0.85, fc = MQ.l,M = = 1,(,/N = 0.3125,;/ = 0.2 
and Ai = fc for both (a) and (b). The regions of stable syn- 
chronized states are indicated by (S). 



each individual oscillator has to be taken into account. 
Namely, for a large relative modulation Xi/k ^ 1 not 
only the individual oscillator frequency is affected but 
also the amplitudes. 

First, we proceed by considering the spatiotemporal 
state of the system (|H)-([n|) with and without external 
parametric excitation. The time evolution of = 100 
self-oscillators from an initially inhomogeneous collective 
state is followed by plotting space-time diagrams. These 
diagrams display a dynamical collective state of the sys- 
tem. In order to do so, we integrate 2N + 1 equations. 
The initial conditions are prepared by adding random 
perturbations to the initial zero state. Figure ^ shows 
the space-time diagram of the amplitude of ith oscilla- 
tor \zi\ on the {i,t) plane for TV = 50 oscillators. The 
oscillator number is shown in abscissa and time t in ordi- 
nate. Figure ([5^) displays the time evolution of an asyn- 
chronous state from the initial time t = to the final 
time t — SOOtt in the absence of external periodic driv- 
ing. A closer look reveals that the elements in the half of 
the ensemble (i = 1, . . . 50) are grouped in small clusters 
of oscillators. Each group consists of elements that are 
phase shifted with respect to the neighboring elements of 
the same group. As the external stimulation turned on 
the state in Fig. ([5^) evolves to a new asynchronous state 
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FIG. 7: Asynchronous state of oscillators near the 1 : 1 reso- 
nant driving region from Fig. ((GJd). (a) Amplitude of total 
extension Xtot = X]i=i^i ^ chain of N = 100 oscilla- 
tors versus external frequency Sle. (b) Two frequency ratios 
versus external frequency are displayed: Q^iv'^^ (dots) and 
Qei^^ot (circles). Um is the average frequency of oscillators in 
a chain, utot in the frequency of total extension, (c) Phase 
= arg(2;i) of TV = 100 active elements as a function of 
time for asynchronous state solution for Sle = 0.85. (d) Ex- 
ample of quasiperiodic solution plotted for xtot versus force 
F. Qe = 0.85 for (c) and (d). 



shown in the panel (b) of the same figure. It presents 
the time-evolution of an inhomogeneous unstable state 
for TV = 50 oscillators and the frequency of excitation 
He = 0.85. It can be easily seen that the distribution of 
groups of oscillators remains the same as in Fig. ([5^) but 
the average self-oscillator frequency is shifted towards 
the stimulation frequency Og. The state is not stable, 
and is sensitive to variations of initial distributions of 
the positions and the velocities of the oscillators. In fact, 
the value of numerically estimated Lyapunov exponent 
is positive in this case. Such regimes are found for the 
CGLE with an additive periodic forcing These "tur- 
bulent synchronized states" arise in the situation when 
the forcing is weak and its frequency is chosen nearly 
the same as a natural oscillation frequency [§]. In this 
case the dynamical state is characterized by the average 
frequency of collective excited oscillations being equal to 
external frequency and its numerical Lyapunov exponent 
being greater than zero. Finally, a perfectly synchronized 
solution is displayed in Fig. Unlike in the previous 

two plots, here after a short transient time the ensemble 
is organized into coherent structure with all the individ- 
ual elements moving in phase and with the frequency fie- 
To produce this plot another type of stimulation is used 
that will be discussed in the following sections. 



We will specify the transitions that occur in the sys- 
tem for various stimulation frequencies by plotting the 
frequency ratio curves versus the stimulation frequency. 
Figure (0 displays two cases for distinct parameter 
choices. Frequency ratio curves for an initially antisyn- 
chronous state and for an initially asynchronous state are 
shown in Fig. ([6]) (a) and (b). Let us first examine the 
changes in dynamical behavior if we chose the parameters 
from the domain that corresponds to an antisynchronous 
oscillatory regime stable for the time-independent sys- 
tem. This can be easily done by defining the bifurcation 
parameter e = 0.1 and by taking the frequency of the in- 
ertial oscillator Hq = 0.85 We plot in Fig. ^) the 
frequency ratio curves fleUJ^^ versus fie- Several distinct 
regions that correspond to different collective dynamical 
states exist. The system remains in the antisynchronous 
state for an interval G [0.2; 1.2]. For higher frequencies 
fie € [1.22; 2.03] the solution is driven into the synchro- 
nized 2 : 1 state. This synchronized state is character- 
ized by all the oscillatory elements moving in-phase and 
with the same frequency. Upon increasing the excitation 
frequency fig > 2.03 the system undergoes a transition 
to an antisynchronous behavior. We examine now the 
transitions in frequency that occur for an initially asyn- 
chronous state. By adjusting the bifurcation parameter 
to e = 0.5 one insures the existence of an asymptoti- 
cally stable asynchronous state of the unperturbed sys- 
tem (Ai = 0) ([4])- ([6]). The frequency of the inertial 
oscillator is chosen to be the same as defined above. We 
present results of numerically calculated frequency ratios 
in Fig. ^p). Two curves are presented: the frequency 
ratio fleUJ^gl between external driving frequency and the 
frequency of the total extension Xtot = J2i=i versus 
fie (squares) and the frequency ratio fleUj:^^ versus fie 
(dotted curve), where ujm is the spatially averaged fre- 
quency. One can observe that transitions from one collec- 
tive state to the other are controlled by changing the fre- 
quency fie- In particular, there exists a 1 : 1 plateau that 
indicates the region of quasiperiodic asynchronous solu- 
tions. There are several transient regimes calculated for 
different values of the external driving fie '■ asynchronous 
states for fie = [0.2; 1.2], the region of coexistence of syn- 
chronous and asynchronous states for fie = [1.2; 1.35] at 
the border with the 2 : 1 resonance. A large set of syn- 
chronized solutions emerges for the 2 : 1 ratio when all 
the oscillatory elements are phase locked to the 2 : 1 fre- 
quency ratio regime. Note, that although for the initially 
antisynchronous state the one-to-one stable entrainment 
regime is not present (see Fig. ([6^)) this regime exist for 
the initially asynchronous state. For a detailed view on 
this regime we refer to Fig. ([7]) which shows a magnifica- 
tion of the 1 : 1 plateau. Numerically calculated nonlin- 
ear amplitude response (see Fig. ([7^)) for a total chain 
extension Xtot and the 1 : 1 plateau in Fig. ((TId) are mag- 
nified for ranges of frequencies fie that lie close to the 
1 : 1 resonance. In Fig. ^) the evolution of oscillator 



8 



0.8 



Ii:3i2:5 111:2 112:3 



U:1; 



0.6 



0.4 



0.2 



0.4 



0.6 

e 



0.8 



FIG. 8: Arnol'd tongues: stability diagram for the parameter 
space defined by amplitude A2 of external driving and fre- 
quency ratio QeLOQ^. The boundaries separating stable limit 
cycles from unstable and quasiperiodic solutions are shown. 
The frequency ratios of entrained limit cycle solutions inside 
each region are indicated. Parameters are: e = —0.2, v — 
0.2,^ = 0.3125,5 = = 1,M = = 0.4. Frequency of 
an unperturbed limit cycle solution ujq = 0.6337. 



phases (pi = arg{zi) is plotted for N = 50. An example of 
quasiperiodic solution for 1 : 1 resonance is presented in 
Fig. JTJi). Here, the oscillation of tension force F is dis- 
played versus the total chain extension xtot- Our results 
suggests, that the final collective behavior depends on the 
choice of initial collective state of the system. Indeed, as 
we have seen in Fig. ^ the same stimulation frequency 
results in two distinct final states. One interesting detail, 
that the stability diagrams indicate the largest set of the 
2 : 1 entrained solutions. The emergence of these solu- 
tions is due to the presence of the multiplicative periodic 
modulation that produces the 2 : 1 harmonic terms in 
the system equations. 



SYSTEM WITH PARAMETRIC MODULATION 
OF FREQUENCY OF LIMIT CYCLE 
OSCILLATORS 

In this section we discuss the synchronization dynamics 
for the system of self-oscillators with a parametric mod- 
ulation of the frequency of self-oscillatory element. We 
modify a set of equations Q-® by setting Ai = and 
introducing a periodic term g[X2,t) = A2Cosf2et with 
the frequency fig and the amplitude A2. The resulting 
equation for the zth oscillator can be written as follows: 
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FIG. 9: Limit cycle solutions for the system (fTOll - (fTT|l . Thick 
line shows the original limit cycle for the autonomous sys- 
tem (|10p - (|lip (A2 — 0). Two limit cycle solutions (thin and 
dotted curves) are displayed for the system pop -pi p with 
X2 = uJ and external frequency = 0.6337 (dotted line), 
and with A2 ~ O.Ilo (thin line). Parameters are the same as 
in Fig. dHl). 



Where g{X2,t) — uj + g{X2, t) is the time-dependent fre- 
quency modulation term. Consequently, the equation for 
the force yields: 



N 



F : 



^F , ^ j^^J d{B\zi\'^z, + ig{X2,t)z. 



NM N ■ 



dt 

{e + v)z,-^llz A, {11) 



z^ = MA2, i) + t)Zi - Sjzipz, + FC 



(10) 



Stability diagram for the synchronized state 

Analogously to the previous section where we first dis- 
cussed the entrainment behavior for the low-dimensional 
system, here we study first the case of = 1 oscillator. 
By choosing e < and the linear oscillator frequency 
JIq below the critical value fie we assert that the sta- 
ble globally synchronized state exist for the autonomous 
system P^ - ([TT|) (A2 — 0). We intergate the set of equa- 
tions pTI)l - (fTT|) with initial conditions picked randomly 
near to the trivial solution. We carry out numerical sta- 
bility for a non-zero parametric modulation g(X2^i\ In 
order to reach the limit of a strong excitation we define 
the amplitude of the modulating term to be comparable 
with the amplitude of a self-oscillating element frequency 
A2/CLI < 1. Figure ([8]) shows Arnol'd tongues calculated 
numerically. Inside the stability regions the frequency 
ratios are identified for each resonance. Examples of the 
limit cycle solutions for the unperturbed system (A2 = 0) 
and for the perturbed system with X2UJ~^ = 0.1 and 
A2i^^^ = 1 are shown in Fig. 
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FIG. 10: Amplitude-frequency response surface. Numeri- 
cally calculated maximum of solution x{t): max|a;| for the 
system (I10|l - (|lll ). The stability boundary of 1 : 1 entrain- 
ment region is shown on the (f2e,A2) plane. Parameters: 
e = -0.2, ly = 0.2, k = 0.4, M = 1, ^ = 0.3125. The frequency 
of an unperturbed limit cycle is uio ~ 0.6337. 



FIG. 11: Frequency ratio curves versus external frequency Qe 
calculated from the Eqs. (|10p and (lllll . Parameters: (e = 0.1) 
for antisynchronous state in (a), (e = 0.5) for asynchronous 
state in (b), Qo = 0.85 and A2 — 0.5lj. 



Amplitude— frequency response 

In Fig. ITUl) the amplitude-response surface is pre- 
sented from numerical calculations of the maximum of 
solution x{t) for the system pHll - dTT]) . The maximum 
is estimated over the entire interval of integration and 
obtained from the zero-averaged solution max|a;| — 
max2007r<t<3007r \x{t) — x\. The borders of stability re- 
gion are displayed on the (Jle, A2) plane. The response is 
smooth inside the one-to-one entrainment stability zone. 
It strongly shows nonlinear features for large amplitude 
of the periodic driving: A2 < w. 



Entrainment for inhomogeneous states 

Now we turn to the case of entrainment behavior for in- 
homogeneous states. We compute frequency ratio curves 
for initially inhomogeneous states of a chain of TV = 100. 
The stimulation frequency is chosen within the interval 
[0.2; 1] in order to include the 1 : 1 resonance and the 
lower ratio resonances. The results are reproduced in 
Fig. pT|) for an initially antisynchronous state (a) and 
for an initially asynchronous state (b). In both cases, one 
can distinguish two well-pronounced extended plateaus 
at the 2 : 3 and the 1 : 1 resonances. The plateaus 
indicate the perfectly entrained solutions with all the os- 
cillatory units moving in phase and with the frequency 
of rig. The remaining solutions that does not belong to 
the entrainment regions are inhomogeneous states. They 
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FIG. 12: Synchronized state of oscillators inside the 1 : 1 
resonant region in Fig. (|llb). External driving frequency 
Q.e = 0.85 for (c) and (d). 



display qualitatively different spacial dynamics. Each 
state is characterized by the individual oscillators form- 
ing macroscopic clusters. The motion of the individual 
units in a cluster is quasiperiodic. One can observe a 
region of multistability for low values of ilg. Stable co- 
existing homogeneous and inhomogeneous states can be 
found inside this region. In the following figure the 
magnification on the one-to-one synchronized state is 
displayed. This state is obtain from an initially ran- 
domly distributed phases of individual elements in the 
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chain. Th amplitude of modulation is A2 = O.Sw and the 
external frequency fie = 0.85. The space-time plot of 
the synchronized state ca be also viewed in Fig. ([St) . Af- 
ter a short transient all the oscillators begin to move in 
phase and with frequency equal to Qe. For the synchro- 
nized state shown in Fig. the frequency of the total 
extension xtot and the frequencies of self-oscillatory ele- 
ments math the stimulation frequency fig- In Fig. (112b .) 
and (b) the amplitude response and the frequency ratio 
are plotted versus fig. The periodic orbit in the space of 
variables F and xtot is also shown. 



DISCUSSION 

In this work we study entrainment behavior in an en- 
semble of identical limit cycle oscillators coupled via geo- 
metrical constraint and via a common force to an inertial 
linear oscillator. The external force due to a mass-spring 
load provides a global coupling: it acts on each element of 
the ensemble. We consider an external parametric modu- 
lation as a time-dependent force acting in addition to the 
applied force due to a mass load The collective dynamics 
of oscillatory elements undergoes several transitions be- 
tween different inhomogeneous states. Such transitions 
are ruled by the frequency and the type of parametric 
modulations introduced in the model. Based on our de- 
scription we find families of stable oscillatory synchro- 
nized solutions corresponding to the resonant external 
stimulation of the system. We discuss an influence of 
two types of parametric excitations: periodic modula- 
tion of the stiffness of the inertial load and modulation 
of the frequency of self-oscillatory element. In particular, 
we demonstrate that for the spatially synchronized state 
the second type of excitation leads to larger stability do- 
mains than the first type of excitation. Furthermore, we 
show that with the appropriate parametric modulation 
introduced an initially inhomogeneous state of the sys- 
tem is driven into the 1 : 1 synchronized state with the 
frequency of collective oscillations equal to the external 
frequency. Our conclusion is that such globally synchro- 
nized behavior can not be achieved from an initially inho- 
mogeneous state by the periodic modulation of the stiff- 
ness coefficient. Instead a new type of dynamical state is 
observed. The emergent spatially inhomogeneous behav- 
ior can be described as a partially synchronized state: 
while every individual element moves with a frequency 
higher than the frequency of the total chain extension, 
the sum of displacements of all the elements oscillate at 
a frequency equal to the external stimulation frequency. 

The amplitude-frequency characteristic shows a typ- 
ical nonlinear response behavior with a fold in the pa- 
rameter space defined by the driving frequency and the 
amplitude of the response. Our numerical results show 
a good agreement with the quasiperiodic approximation 
inside the one-to-one stable entrainment zone. 



The phenomena discussed here have relevance and ap- 
plications in nature and laboratories. In particular, the 
autonomous model and its time-dependent analogue pro- 
vide a phenomenological approach to model oscillatory 
dynamics of an active muscle fiber. Experimental in 
vitro observations on skinned muscle fiber show evidence 
of different oscillatory regimes generated by spontaneous 
contractions of muscle elementary units, sarcomeres [2ll |. 
Depending on the conditions of the experiment, oscil- 
latory contractions can undergo a synchronized activ- 
ity [2l| or a spatially inhomo gene ous (asynchronous or 
antisynchronous) activity [25I. l2d|. In experiments the 
signalling between adjacent sarcomeres is chemically reg- 
ulated. The state of activation of one sarcomere is af- 
fected by the state of the adjacent sarcomere 2l|. Our 



current model does not consider the effect of coupling 
between adjacent sites in the chain of elements. In the 
future we plan to include local coupling between ele- 
ments and to consider its effect on spatiotemporal dy- 
namics. Detailed study of the corresponding regimes of 
spatiotemporal behavior is left to future work. 

Experimental work and results from theoretical analy- 
sis of our phenomenological model promise new interest- 
ing insights into collective behavior of oscillations in iso- 
lated fibrillar muscle. Periodic modulation of the model 
parameters considered here might serve as an appropri- 
ate setting for modeling active muscle contractions under 
a specially designed external mechanical and chemical 
stimulation. 



APPENDIX 

In the appendix, we discuss a quasiperiodic approxima- 
tion 



27| used for derivation of the nonlinear amplitude- 



frequency dependence. The derivation is carried out for 
the first type of parametric modulation discussed in the 
paper. It can be seen from Fig. ([2]) that the original 
unperturbed limit cycle is close to the ellipse centered at 
zero. Solutions for a non-zero small perturbation Ai <^ fc 
are nearly close to the disturbed elliptical form with their 
centers shifted away from zero. These observations lead 
us to seek the solutions expansion in form of a finite series 
of harmonic terms, where the largest contribution terms 
represent solution for an ellipse. The system dU-dS]) can 
be rewritten in a real form after introducing x = Re(z) 
and y = Im(z) variables with x a displacement of the 
oscillator and y a time-dependent variable nonlinearly 
coupled to x: 



M/^x = -{l + Mv/^)x-K{\i,t)/^x 

-Ljy + ex - B{x'^ + y'^)x, (12a) 
y = ujx + ey- B{x^ +y^)y, (12b) 
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For the non-autonomous system ((12)) a limit cycle solu- 
tion can be written in the quasiperiodic form |28] : 



X = ai cos fit + bi sin Qt + uq 
y = b2 sin fit + 02 cos fit + ag 



(13a) 
(13b) 



Where fl is the frequency of a new limit cycle for the 
time-dependent system ([T^. The amplitude coefficients 
in the expansion (fT3)) depend on the size of the orig- 
inal unperturbed limit cycle solution and are ordered 
according to the powers of e: ai,b2,ag are of order 
Oie^/"^), 6i,a2 areof orderO(e) andao = 0(e3/2). j^^j^g 
driving frequency fl,, is near to the frequency loq of the 
limit cycle solution for the unperturbed system (Ai — 0) 
in (III)-® a reasonable approximation for the 1 : 1 lock- 
in solution can be made by retaining only terms to the 
accuracy of O(e^) in Eqs. (fT3|) . We would not consider 
the terms of order equal to O(e^) in the approximation 
presented here. Coefficients of the expansion (|13p can be 
expressed via coefficients in ^ as follows: 



a-i = , oi ^ , ao = a;, (14a) 



0-2 



62 = 



i(d_i - di) 



ig = y. (14b) 



We introduce a small parameter A = f2 — fig = O(e^) to 
be a frequency mismatch between the external frequency 
and the resulting limit cycle frequency. Let us define a 
linear frequency: 

r!o = w/(l + A'WC)i/2. 

Consequently, one finds that the frequency ratio fleUJ^^ 
can be expressed as follows: 

fl^u^'^ w flfl-'^{l ~ A). 

In the next step, we use the quasiperiodic approxi- 
mation (fT5)) to derive a final form of the amplitude- 
frequency response. We assume that the following ex- 
pansion holds true: A: = /co + efci, with fco = Mfl^. For 
the sake of simplicity, we introduce the rescaling of the 
parameters: 



Upon substituting the quasiperiodic form ()13|) into 
Eq. ([5]) and collecting various terms according to the or- 
der of harmonics we obtain a set of nonlinear equations: 



(15a) 
(15b) 



md'o = — (1 + mv)di) + (e — k)aQ 
—ujag — Aiai/2 — 51 ao 
~aig2 - bigs, 
a's = + - aggi ~ 0232 

-^233, 

d'l = —(1 + 'mv)di + fP'rhai 

— (1 + mv)flbi + (e — k)ai — uja2 

-Aiflo - fligi - 2aog2 

-ai54/4 - 6155/2, (15c) 

0,2 — — f2&2 + ^0L2 + — 2ag52 

-0231 - 02.94/4 - 62.95/2, 

62 = f^a2 + £62 + w6i - 2ag53 - 6251 

-02.95/2 + 6234/4, (15d) 
rhbi = — (1 + mi^)6i + ?7if2^6i 

+ (1 + mv)flai + (e — k)bi — ujb2 

-2ao53 - 0135/2 - 6151 

-6154/4, (15e) 

where 51, 52; 53 are nonlinear functions of the coefficients 
in m 



5i 




2 0? 


a, 

+ y 


52 


= B{aoai 


+ 0902), 




53 


= B{aobi 


+ 0362), 




54 


= B{al^ 


bj + al- 




55 


= B{aibi 


+ 0262). 





6| 
2 



2 



(16) 

We would not consider the terms of order 0(e^/^) in the 
approximation presented here. Let us make the assump- 
tion that for the entrained limit cycle solution ([TH]) the 
coefficients in the expansion (jl3p are constant or vary 
slowly with time. Then one can neglect derivatives from 
the equations for the coefficients (IT5t . By neglecting 
terms of order 0{e^^'^) one can express the coefficients 
a2,ag, 61 and 62: 



'2uj 



02 



61 



rri(flQ 



-ai. 



(17a) 
(17b) 

(17c) 



Let us use the substitution p ^ a\. Our final goal is to 
obtain a family of amplitude-frequency curves: 

T^^^{fl,) : n,u;^^ pifin^'). 

After substituting the expressions (fT7)) for the coefficients 
into Eqs. (fT5|) and eliminating terms of the order 0(e^/^) 
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we write down the approximate amplitude-frequency re- lation: 
I 




Finally, what remains to be done is to determine a set 
of parameters for calculating the relation Let us 

define the parameters in order to compare with the nu- 
merical results: e = — 0.2,j/ = 0.2, Af = l,uj — 1,C = 
0.3125, = 0.7861, Ai = fco,ao = O.OlAj^^ We have 
used A — 0.07 rio^^^^, fci — 0.4^ in order to obtain a 
reasonable fitting to the numerical amplitude-frequency 
response in Fig. 
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